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Abstract 


Chondrules are millimeter-sized spherules that dominate primitive meteorites (chondrites) originating from the asteroid belt. 
The incorporation of chondrules into asteroidal bodies must be an important step in planet formation, but the mechanism is not 
understood. We show that the main growth of asteroids can result from gas-drag-assisted accretion of chondrules. The largest 
planetesimals of a population with a characteristic radius of 100 km undergo run-away accretion of chondrules within ~3 Myr, 
forming planetary embryos up to Mars sizes along with smaller asteroids whose size distribution matches that of main belt asteroids. 
The aerodynamical accretion leads to size-sorting of chondrules consistent with chondrites. Accretion of mm-sized chondrules 
and ice particles drives the growth of planetesimals beyond the ice line as well, but the growth time increases above the disk life 
time outside of 25 AU. The contribution of direct planetesimal accretion to the growth of both asteroids and Kuiper belt objects is 
minor. In contrast, planetesimal accretion and chondrule accretion play more equal roles for the formation of Moon-sized embryos 
in the terrestrial planet formation region. These embryos are isolated from each other and accrete planetesimals only at a low rate. 
However, the continued accretion of chondrules destabilizes the oligarchic configuration and leads to the formation of Mars-sized 
embryos and terrestrial planets by a combination of direct chondrule accretion and giant impacts. 


Introduction 

The formation of planetesimals and planetary embryos is an 
important step towards the assembly of planetary systems ( 1 , 2 ). 
The asteroid belt contains planetesimals left over from the for¬ 
mation of the solar system and thus provides a record of the 
early stages of planet formation. Chondrite meteorites are 
fragments of asteroids that did not undergo large-scale melt¬ 
ing and differentiation. Their constituent particles allow us to 
study the first stage of planet formation, where micrometer¬ 
sized dust grains grew to asteroid-scale planetesimals. The 
dominant mass fraction of most chondrite meteorites are chon¬ 
drules, millimeter-sized glassy spherules formed by transient 
heating events in the protoplanetary disk (5). The small ma¬ 
trix grains that fill the space between these chondrules likely 
entered the meteorite parent bodies as accretionary rims on the 
chondrules ( 4 ). 

Uranium-corrected Pb-Pb dates of chondrules show crys¬ 
tallisation ages ranging from the earliest epoch of the solar 
protoplanetary disk - defined by the condensation of calcium- 
aluminium-rich inclusions (CAIs) 4567.30 ± 0.16 Myr ago - 
to approximately 3 Myr later (5). Moreover, chondrules from 
individual chondrites show variability in their 54 Cr/ 52 Cr and 
50 Ti/ 1 ' Ti ratios (5, 6 ) that track genetic relationships between 
early-formed solids and their respective reservoirs. Collec¬ 
tively, these observations indicate that chondrules from individ¬ 
ual chondrite groups formed in different regions of the proto¬ 
planetary disk and were subsequently transported to the accre¬ 


tion regions of their respective parent bodies. Thus, the region 
of asteroid formation must have been dominated by chondrules 
over time-scales comparable to the observed lifetimes of proto¬ 
planetary disks around young stars (7). 

Meteorites provide a direct view of the primitive, chon- 
dritic material from which the planetesimals in the asteroid 
belt formed. In contrast, planetesimals originally located in 
the formation region of terrestrial planets accreted into plan¬ 
etary bodies and hence little evidence exists with respect to 
the material from which those planetesimals formed. However, 
the abundance of chondrules in the asteroid belt suggests that 
chondrules could have been widespread in the terrestrial planet 
formation region as well. Indeed, some chondrules record 
54 Cr/ 52 Cr and 50 Ti/ 47 Ti ratios indicating formation in the ac¬ 
cretion regions of Earth and Mars (5, 6 ). Hence, understanding 
the role of chondrules in the formation of asteroids can provide 
critical insights into how the terrestrial planets formed closer to 
the Sun. 

The mechanism by which large amounts of chondrules were 
incorporated into asteroids is poorly understood. Particles of 
chondrule sizes can concentrate near the smallest scales of the 
turbulent gas ((S’). The most extreme concentration events have 
been proposed to lead to asteroid formation by gravitational 
contraction of such intermittent chondrule clouds ( 9 ). How¬ 
ever, it is still debated whether such high concentrations actu¬ 
ally occur ( 10 ). The streaming instability ( 11 ) is an alternative 
mechanism that concentrates larger particles ( 12 , 13 ), of typical 
sizes from 10 cm to 1 m when applied to the asteroid belt, due 


1 





UH 


Figure 1: The maximum particle density versus time for 
streaming instability simulations without self-gravity at resolu¬ 
tions 64 3 (black line), 128 3 (yellow line) and 256 3 (blue line). 
The particle density is measured in units of the mid-plane gas 
density p g and the time in units of the orbital period T or b. The 
first 30 orbits are run with a solids-to-gas ratio of Z = 0.01, 
with only modest overdensities seen in the particle density. 
Half of the gas is then removed over the next 10 orbits, trig¬ 
gering strong particle concentration, of up to 10,000 times the 
local gas density at the highest resolution. Doubling the res¬ 
olution leads to a more than a quadrupling of the maximum 
particle density. 

to an aerodynamical effect where particles pile up in dense fil¬ 
aments, which can reach densities of more than 1,000 times the 
local gas density (14). The streaming instability nevertheless 
fails to concentrate chondrule-sized particles whose motion is 
strongly damped by gas drag (15). 

Results 

Here we report the results of a model where asteroids and 
planetary embryos grow by accretion of chondrules onto plan- 
etesimal seeds formed by the streaming instability. Particles 
of large enough size to concentrate by streaming instabilities 
could have been present in the early stages of the protoplan¬ 
etary disk when planetesimals formed. Such particles include 
macrochondrules (16), chondrule aggregates (17) and ice-rich 
pebbles. All of these drift rapidly towards the Sun because of 
gas drag. Therefore we assume that cm-sized particles were 
only present in the earliest stages of the protoplanetary disk, 
whereas smaller chondrules existed during most of the disk’s 
life-time. Chondrule-sized particles may have been able to re¬ 
main in the asteroid belt due to stellar outflows (18) or advec- 
tion with the outwards moving gas in the mid-plane layer of 


Figure 2: The maximum particle density versus the length scale 
(measured in scale-heights H). We have taken spheres with di¬ 
ameters from one grid cell up to the largest scale of the sim¬ 
ulation and noted the maximum value of the density over all 
simulation snapshots. The results are relatively converged on 
each scale. The increase in maximum particle density with in¬ 
creasing resolution is an effect of resolving ever smaller-scale 
filaments. Blue dotted lines mark the Roche density at three 
values of the gas column density at 2.5 AU compared to the 
Minimum Mass Solar Nebula (MMSN). The red dotted line in¬ 
dicates the characteristic length scale of the streaming insta¬ 
bility. The black line shows a power law of slope —2, which 
shows that the maximum density follows approximately the in¬ 
verse square of the length scale. 

sedimented particles (15). We discuss the conditions for plan- 
etesimal formation in the asteroid belt, particularly the forma¬ 
tion of dm-sized particles, further in the Supplementary Text. 

Planetesimal formation simulations 

We use high-resolution models of planetesimal formation 
through streaming instabilities to set the initial conditions for 
our model. The highest resolution reached prior to this work in 
such simulations is 128 3 grid cells with 2.4 million superparti¬ 
cles representing the solids (14). We here report on simulations 
with up to 512 3 grid cells and 153.6 million superparticles rep¬ 
resenting the solids, which is 64 times more resolution elements 
than in previous simulations. Details of the simulation method 
are given in Materials and Methods. 

We first perform a convergence test of streaming instability 
simulations with no self-gravity between the particles. Figure]!] 
shows the maximum particle density versus time in these sim¬ 
ulations. We run the initial 30 orbits with the full gas density. 
Here the maximum particle density is very moderate and no 
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strong concentrations occur. Between 30 and 40 orbits we re¬ 
move 50% of the gas, effectively increasing the solids-to-gas 
ratio from Z = 0.01 to Z = 0.02. This increase in metallic - 
ity triggers strong particle concentration through the streaming 
instability, of up to 10,000 times the local gas density. The 
maximum particle density increases approximately quadrati- 
cally with the inverse size of the resolution element Sx. The 
maximum particle density on different scales is shown in Fig¬ 
ure^ Here we have measured the maximum density in spheres 
of diameters from Sx up to the full size of the simulation box. 
The results are well converged from 64 3 to 256 3 at the scales 
that are present at all resolutions. The smaller scales made ac¬ 
cessible at higher resolution yield increasingly higher densities. 
This is an effect of the filamentary structure of the overdense 
particle filaments forming by the streaming instability; higher 
resolution allows us to resolve finer structure and hence higher 
densities. Smaller planetesimals can thus form as the resolution 
is increased and increasingly smaller-scale filaments cross the 
Roche density. 

In Figure[3]we show size distributions of planetesimals form¬ 
ing in streaming instability simulations that include the self¬ 
gravity of the particles. The planetesimal size distribution is 
dominated in number by small planetesimals, but most of the 
mass is in the few largest bodies. We find that progressively 
smaller planetesimals form, alongside the large planetesimals, 
at higher resolution. The largest planetesimals, which con¬ 
tain the dominant mass of the population, decrease to approx¬ 
imately 100 km in radius for the column density of the Mini¬ 
mum Mass Solar Nebula at the location of the asteroid belt (79). 
The size distribution of asteroids above 60 km in radius has 
been shown to retain its primordial shape, since the depletion 
of the asteroid belt happens in a size-independent fashion for 
those large bodies (20,21). The differential size distribution re¬ 
sulting from streaming instabilities disagrees with the observed 
size distribution of the asteroids (dashed line in Figure[3]). 
Chondrule accretion on asteroids 

We proceed to consider the subsequent evolution of the size 
distribution as the newly formed planetesimals accrete chon- 
drules present in the gaseous surroundings. We assume that 
planetesimals are born within an ocean of small chondrules 
that do not directly participate in planetesimal formation. Ac¬ 
cretion of these chondrules by the planetesimals over the fol¬ 
lowing millions of years can lead to significant mass growth. 
We use initial size distributions inspired by the streaming in¬ 
stability simulations. However, chondrule accretion is not lim¬ 
ited to this model; it would be efficient in the context of any 
planetesimal formation mechanism that produces planetesimals 
with characteristic sizes of 100 km (9,22). 

The gas in the protoplanetary disk is slightly pressure- 
supported in the radial direction and hence moves at a sub- 
Keplerian speed, with Av marking the difference between the 
Keplerian speed and the actual gas speed. The speed difference 
is approximately 53 m/s in the optically thin Minimum Mass 
Solar Nebula model (79). The Bondi radius of a planetesimal 
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Figure 3: Birth size distribution of planetesimals forming by 
the streaming instability in 25-cm-sized particles. The differ¬ 
ential number of planetesimals (per 10 22 g) is calculated with 
respect to the nearest size neighbours in the simulation. Yellow, 
blue and red circles indicate individual planetesimals forming 
in computer simulations at 128 3 , 256 3 and 512 3 grid cells. The 
size distribution of the highest-resolution simulation is fitted 
well with a power law dN/dM oc M~ qM (red line; the fit is 
based on the cumulative size distribution shown in Figure[4] but 
does not include the exponential tapering). Simulations with 
lower values of the particle column density S p yield succes¬ 
sively smaller radii for the largest planetesimals, down to below 
100 km in radius for a column density similar to the Minimum 
Mass Solar Nebula model (27 p = 4.3gem -2 at r = 2.5 AU). 
Binary planetesimals (marked with circles) appear only in the 
highest-resolution simulation, as binary survival requires suf¬ 
ficient resolution of the Hill radius. The differential size dis¬ 
tribution of the asteroid belt (dashed line) shows characteristic 
bumps at 60 km and 200 km. The planetesimal birth sizes from 
the simulations are clearly not in agreement with main belt as¬ 
teroids. 


with radius R and internal density p,. 

Rb _ ( R \ 2 ( Av \- 2 / p. 

R \ 100 km/ \ 53 m/s/ \3.5g/cm 3 

measures the gravitational deflection radius of the planetesimal. 
Chondrules embedded in the sub-Keplerian gas flow are ac¬ 
creted from the Bondi radius (23,24) when their friction time if 
is comparable to the time to drift azimuthally across the Bondi 
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Figure 4: Cumulative birth size distribution of planetesimals forming by the streaming instability in our highest- 
resolution simulation (512 3 ). The cumulative size distribution is fitted with an exponentially tapered power law A’> oc 
M -0 - 6 exp[—(M/M exp )( 4 / 3 )] (red line), with exponential tapering at M exp = 1.2 x 10 23 g (dotted line). Shallower or steeper 
power laws yield poorer fits to the populations of small and large planetesimals, respectively. We choose to fit the cumula¬ 
tive size distribution rather than the differential size distribution to avoid the noise inherently present in the latter. The fit can 
be translated to dN/dM oc M -16 exp[— (M/M exp ) 4 / 3 ] (the power law part of this fit is indicated in Figure [ 3 J as well as to 
dN/dR oc R~ 2 8 exp[—(f?/f? exp ) 4 ] in the differential size distribution per unit radius. 


radius, t\>, = Rb/Av. An additional turbulent component to 
the gas motion can be ignored since this is expected to be much 
slower than the sub-Keplerian speed in the asteroid formation 
region. Peak accretion rates are obtained for ff/fe in the range 
from 0.5 to 10 (24). This happens at the orbital distance of the 
asteroid belt for particle sizes in the interval 

«=[0.08,1.6]mm( Is y (J/Nj /»„. (2) 

Here / gas represents the ratio of the actual gas column density 
to that of the Minimum Mass Solar Nebula (19). Chondrule 
accretion happens at a rate proportional to /i| 5 (or M oc R 6 ). 
This run-away accretion can drive a very steep differential size 
distribution similar to what is observed for large asteroids in the 
asteroid belt. 

We solve numerically for the temporal evolution of the sizes 
and orbits of planetesimals accreting chondrules and collid¬ 
ing with each other in an annulus extending from 2.4 to 2.6 
AU (see Materials and Methods for details). The planetes¬ 
imals have initial radii from 10 to 150 km, distributed in 
size along an exponentially tapered power law based on the 
cumulative size distribution of planetesimals forming in the 
streaming instability simulations (see Figure |4j), dN/dR oc 
R~ 2 S exp[—(f?/f? exp ) 4 ], with exponential tapering at R exp = 
100 km and a total mass of 0.04 Earth masses. Chondrules have 
diameters from 0.02 mm to 1.6 mm, in broad agreement with 
observed chondrule sizes (3). The growth of particles and the 
formation of chondrule precursors is discussed in the Supple¬ 
mentary Material. We place 0.2 Earth masses of chondrules 


initially in the annulus while another 0.2 Earth masses is cre¬ 
ated continuously over a time span of 3 Myr. The chondrules 
are given a Gaussian density distribution around the mid-plane, 
with the scale-height set through a nominal value of the turbu¬ 
lent viscosity of a = 10 . The gas in the protoplanetary disk 

is depleted exponentially on a time-scale of 3 Myr (7). 

We show in the left panel of Figure|5]the size distribution af¬ 
ter 5 Myr of accreting chondrules. The steep drop in the initial 
number of planetesimals larger than 100 km in radius has be¬ 
come much shallower because of chondrule accretion, forming 
a characteristic bump at 60-70 km radius. The size distribution 
of smaller asteroids in today’s belt was likely filled later with 
fragments from collisional grinding (20). The steep size distri¬ 
bution of larger asteroids ends at around 200 km as even larger 
asteroids would need to accrete chondrules larger than millime¬ 
ter in size (see equation 2). At this point the largest asteroids 
enter a more democratic growth phase where the accretional 
cross sections are reduced by gas friction within the Bondi ra¬ 
dius (23,24). The agreement with the observed size distribution 
of asteroids is excellent in the nominal model, except in the 
range around Ceres mass where it predicts too few objects by 
a factor of approximately 2-3. Changing the model parameters 
(using larger chondrules or smaller planetesimal birth sizes) 
yields poorer agreement with observed asteroid sizes. This im¬ 
plies that the observed size distribution of asteroids is directly 
determined by the size distribution of the chondrules and the 
birth sizes of the planetesimals. 

The right panel of Figure [5] shows that the major growth 
phase of asteroids with final radii ranging from 200 km to 500 
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Figure 5: The size distribution of asteroids and embryos after accreting chondrules for 5 Myr (left panel) and selected masses 
and inclinations as a function of time (right panel). The nominal model (black line) matches closely the steep size distribution of 
main belt asteroids (red line, representing the current asteroid belt multiplied by a depletion factor of 540) from 60 km to 200 km 
in radius. The size distribution becomes shallower above 200 km; this is also seen in the asteroid belt, although the simulations 
underproduce Ceres-sized planetesimals by a factor of approximately 2-3. A simulation with no chondrules (blue line) produces 
no asteroids larger than 300 km. Inclusion of chondrules up to cm sizes (pink line in insert) gives a much too low production 
of Ceres-sized asteroids, while setting the exponential cut-off radius of planetesimals to 50 km (green line) leads to a poorer 
match to the bump at 60 km. The right plot shows that the formation of the first embryos after 2.5 Myr quenches chondrule 
accretion on the smaller asteroids by exciting their inclinations i (right axis). The dotted lines indicate the mass contribution 
from planetesimal-planetesimal collisions. Asteroids and embryos larger than 200 km in radius owe at least 2/3 of their mass to 
chondrule accretion. 


km occurs after 2.5 Myr. Beyond this time the largest planetesi¬ 
mals in the population grow to become planetary embryos with 
sizes similar to the Moon and Mars. These growing embryos 
excite the inclinations of the smaller asteroids to * ~ 0.01, 
which disconnects the asteroids from the chondrules in the mid¬ 
plane layer, quenching their accretion of further chondrules. 
The beginning of embryo formation terminates the accretion 
of asteroids after 3 Myr, thus defining the final sizes of the as¬ 
teroid belt population. The absence of such planetary embryos 
in today’s asteroid belt may reflect a later depletion by grav¬ 
itational perturbations from Jupiter (25-27). Depletion of the 
asteroid belt is discussed further in the Supplementary Text. 
An important implication of chondrule accretion is that accre¬ 
tion of other planetesimals contributes only a minor fraction to 
the final masses of large asteroids and embryos (right panel of 
Figure |5j. 

Size sorting of chondrules 

Chondrules in chondrites appear strongly size-sorted (8,28), 
with the average chondrule diameter varying among the ordi¬ 
nary chondrites (28) from 0.3 mm (H chondrites) up to 0.5 
mm (L and LL chondrites). Carbonaceous chondrites exhibit 


a larger range in chondrule sizes, from 0.1 to 1 mm. In Fig¬ 
ure [6] we show that the mean diameter of accreted chondrules 
in our model lies in a range similar to that observed for chon¬ 
drites. The size distribution of accreted chondrules is very nar¬ 
row, also in agreement with the observed size distribution in 
the ordinary chondrites (28). The chondrule accretion process, 
which leads to aerodynamical sorting of chondrules, may thus 
represent the underlying physical mechanism for the size sort¬ 
ing of chondrules observed in chondrite meteorites. Although 
aerodynamical sorting has been previously proposed to arise 
from the gas flow around asteroids (29), our simulations show 
that accretion from the full Bondi radius is much more efficient 
in achieving aerodynamical sorting of chondrules. 

Terrestrial planet formation with chondrules 

Our identification of chondrule accretion as driving planetes- 
imal growth in the asteroid belt implies that chondrules could 
play an important role for terrestrial planet formation as well. 
To test this hypothesis we have performed a numerical inte¬ 
gration of the evolution of planetesimal orbits and sizes at 1 
AU (Figure [7]l. The left panel shows the size distribution of the 
planetesimals at 0, 1,3 and 5 Myr. Growth to super-Ceres-sized 
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Figure 6: Mean chondrule sizes ((d), upper panel) as a func¬ 
tion of layer radius R and size-distribution width (a 0 , lower 
panel) as a function of mean chondrule size. Yellow lines in 
the upper panel indicate the chondrule size evolution in indi¬ 
vidual asteroids and embryos, while red lines indicate mean 
accreted chondrule sizes at different times. The accreted chon¬ 
drule size increases approximately linearly with planetesimal 
size at t = 0 Myr. Asteroids stirred by the growing embryos 
over the next 2 Myr accrete increasingly larger chondrules, as 
asteroid velocities align with the sub-Keplerian chondrule flow 
at aphelion. Finally, asteroids accrete surface layers of mainly 
very small chondrules, down to below 0.1 mm in diameter, at 
late times when their large inclinations decouple the asteroid 
orbits from the large chondrules in the mid-plane. The width 
of the chondrule size distribution in the lower panel is given in 
terms of a 0 , the base-2-logarithmic half-width of the cumula¬ 
tive mass distribution of chondrules = 1 implies that 2/3 
of the chondrules have diameters within a factor 2 1 = 2 from 
the mean). Dots indicate chondrule layers inside asteroids and 
embryos. Chondrules are aerodynamic ally sorted by the accre¬ 
tion process, to values of a 0 in agreement with the chondrules 
found in ordinary chondrites (hashed region). Unfiltered accre¬ 
tion from the background size distribution of chondrules (blue 
dot: size distribution of un-sedimented particles; yellow dot: 
size distribution in the mid-plane) yields specific pairs of (d) 
and O-0 that are not consistent with constraints from ordinary 
chondrites. 


planetesimals is rapid and happens within 1 Myr. This growth 
is driven mainly by planetesimal accretion, in stark contrast to 


the situation at 2.5 AU. This is seen in the right panel of Figure 
[TJwhere the mass of the largest body is shown as a function of 
time (full line) together with the mass contribution from chon¬ 
drules (dashed line) and from planetesimals (dash-dotted line). 
However, chondrule accretion becomes dominant after 1.5 Myr 
and drives the further growth up to Mars-mass embryos after 4 
Myr. The largest body experiences a giant impact after 4 Myr, 
after which it continues to grow by chondrule accretion towards 
0.9 Earth masses. 

Chondrule accretion can thus promote the growth of the 
largest embryos from Moon masses towards Mars masses and 
finally Earth masses. The dominance of planetesimal accre¬ 
tion in the initial growth towards Moon masses occurs because 
chondrules couple tightly to the gas at 1 AU. Here the gas 
density is more than a factor 10 higher than at 2.5 AU. This 
prevents sedimentation, such that all chondrule sizes are well 
mixed with the gas, and it reduces the cross section for ac¬ 
creting chondrules since tightly coupled particles can not be 
accreted from the full Bondi radius. This situation changes as 
gas dissipates exponentially over 3 Myr, increasing the friction 
times and thus the accretion efficiency of the chondrules. Fur¬ 
thermore, the increasingly large embryos obtain high accretion 
radii for chondrules and hence high chondrule accretion rates, 
despite the relatively low degree of sedimentation of chondrules 
present at 1 AU orbits. 

Chondrules also play a critical role in terrestrial planet for¬ 
mation in a different way, namely by breaking the isolation 
mass configuration. Growth by pure planetesimal accretion 
tends to end in oligarchic growth where the largest embryos are 
isolated from each other by approximately 10 Hill radii (30). 
We have implemented the effect of this isolation tendency by 
identifying the group of isolated bodies as those that can fit 
their combined reach of 10 Hill radii into the annulus of 0.2 
AU in diameter (31, 32). Isolated bodies are not allowed to 
accrete each other and only affect each other dynamically via 
distant viscous stirring in the eccentricity. Inclination pertur¬ 
bations, as well as dynamical friction in the eccentricity, are 
ignored between isolated bodies. 

We mark the isolation mass of approximately 0.01 Me in the 
right plot of Figure [ 7 ] (blue line). A simulation performed with 
no chondrule accretion (red line) shows the expected growth 
of the largest embryos to just below the isolation mass. The 
growth curve for pure planetesimal accretion follows closely 
the growth curve for the full simulation until approximately 
10 -4 Earth masses. At this point chondrule accretion becomes 
significant and the two curves diverge. Importantly, chondrule 
accretion can destabilise a set of isolated bodies, as their con¬ 
tinued growth by chondrule accretion drives the least massive 
of the isolated bodies out of isolation. Hence the giant impacts 
experienced by the largest embryo between 1.5 and 4 Myr are 
all driven by chondrule accretion, as these impacts happen be¬ 
yond the isolation mass. 

The importance of chondrule accretion increases if chon¬ 
drules in the terrestrial planet formation region are larger. 
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Figure 7: Growth of embryos and terrestrial planets at 1 AU. The left panel shows the size distribution at four different times and 
the right panel shows the mass of the most massive body in the simulation as a function of time. The growth up to 1000-km-sized 
embryos is mainly driven by planetesimal accretion, since chondrule-sized pebbles are tightly coupled to the gas and hence hard 
to accrete. However, chondrule accretion gradually comes to dominate the accretion as the embryos grow. The largest body 
reaches Mars size after 3 Myr, with more than 90% contribution from chondrule accretion. A giant impact occurs just before 
4 Myr, where the largest body accretes the third largest body in the population. The continued accretion of chondrules leads to 
the formation of an Earth-mass body after 5 Myr. A simulation with no chondrules evolves very differently: a large number of 
embryos form with masses just below the isolation mass of M- iso « 0.01-Me. 


namely up to 1 cm (macrochondrules). We show the results 
of expanding the size distribution of chondrules up to 1 cm in 
Figure [8] Here the initial growth has equal contribution from 
chondrules and planetesimals. The growth towards embryos 
and terrestrial planets is much more rapid than in Figure [ 7 ] a 
planet with the mass of the Earth already forms after 1.5 Myr. 

The final size of embryos and planets forming by combined 
planetesimal accretion and chondrule accretion depends on the 
exact timing of disk dissipation. If the inner disk already dis¬ 
sipates after 3 Myr, then the terrestrial planet formation region 
will be dominated by a number of Mars-sized embryos. These 
bodies can go on to partake in a traditional terrestrial planet 
formation scenario of gradual orbital perturbations and giant 
impacts over the next 30-100 Myr (55). In contrast, later disk 
dissipation or the presence of large chondrules in the terrestrial 
planet formation region allows for the direct formation of plan¬ 
etary bodies within a few Myr. This is far more rapid than tradi¬ 
tional terrestrial planet formation scenarios that do not consider 
chondrule accretion. 

Pebble accretion in the Kuiper belt 

The size distribution of trans-Neptunian objects has a charac¬ 
teristic bump around 50 km in radius (54), similar to main belt 
asteroids, followed by a steep decline towards larger sizes (20). 
We therefore explore here whether accretion of chondrules and 
icy pebbles could be responsible for the observed size distribu¬ 


tion of Kuiper belt objects as well. 

The asteroid belt displays a value of q « 4.5 in the size 
distribution dN/dR ex R~ q of asteroids larger than 60 km in 
radius, while the hot population of the classical Kuiper belt and 
the Neptune trojans have a much higher value of q « 5 ... 6 
for large planetesimals (55). The cold population of the clas¬ 
sical Kuiper belt has an even steeper decline, with q ^ 8 ... 9. 
The planetesimals in the cold population are believed to have 
formed in situ in the outermost regions of the solar system (56), 
while the hot population and scattered disk objects formed 
outside the birth location of Neptune and were subsequently 
pushed outwards by Neptune’s migration to its current orbit 
(57, 38). The planetesimals which formed between the orbits 
of Jupiter and Neptune may today be represented by the C-type 
asteroids, scattered to the asteroid belt during Jupiter’s migra¬ 
tion (27). In order to probe planetesimal growth by chondrule 
accretion (or icy particles of a similar size range) beyond the 
asteroid belt, we therefore consider the two characteristic dis¬ 
tances 10 AU and 25 AU. 

Planetesimal growth at 10 AU 

The region around 10 AU represents the conditions where 
the gas giants formed. The formation of Jupiter and Saturn 
must have had an enormous effect on the planetesimals in that 
region, with some ejected to the Oort cloud and others poten¬ 
tially mixed into the outer asteroid belt during Jupiter’s and 
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Figure 8: Growth of embryos and terrestrial planets at 1 AU, with chondrule sizes up to 1 cm. The increase in chondrule size 
compared to the previous figure enhances the chondrule accretion rate substantially. The initial growth to 1000-km-sized embryos 
now has approximately equal contribution from planetesimal accretion and chondrule accretion. An Earth-mass terrestrial planet 
forms already after 2 Myr, driven by a combination of chondrule accretion and giant impacts. 


Saturn’s migration (27). To test how the size distribution of 
planetesimals evolves at 10 AU, we set up an annulus of width 
0.2 AU. The planetesimals are again given sizes from 10 to 
150 km in radius, with a steep exponential tapering above 100 
km. The internal density of the planetesimals is set to 2 g/cm 3 , 
similar to that of Ceres, a typical representative of the icy aster¬ 
oids in the outer part of the main belt. The pebbles have radii 
from 0.01 to 1 mm. These pebbles could represent a mixture 
of chondrules transported outwards from the inner solar sys¬ 
tem and icy pebbles formed beyond of the ice line. Chondrules 
may also have formed in situ in the outer parts of the proto¬ 
planetary disk. In fact, some chondrules present in carbona¬ 
ceous chondrites have 50 Ti/ 4 ' Ti ratios comparable to those of 
the most pristine chondrites, namely the Cl carbonaceous chon¬ 
drites (6). This suggests that at least a fraction of chondrules in 
carbonaceous chondrites formed from pristine, thermally un¬ 
processed precursor material typical of the accretion regions of 
Cl chondrites. Given that these chondrites are believed to have 
accreted beyond the ice line (27), chondrule formation appears 
to not have been limited to the asteroidal belt, but also occurred 
in the outer solar system. 

The results at 10 AU are shown in Figure [9] Planetesimals 
of Ceres size grow within approximately 2 Myr, followed by 
the usual phase of embryo growth also seen in Figure [5] The 
largest embryo reaches a mass comparable to the Earth’s be¬ 
fore 3 Myr. The mid-plane density of pebbles is much lower at 
10 AU than at 2.5 AU, but this is more than counter-acted by 
increased sedimentation of particles in the dilute gas. The size 
distribution is very steep, much steeper than in the asteroid belt. 


and resembles better the steepness of the size distribution of 
the trans-Neptunian populations (35). Finally, it is interesting 
to note that at these orbital distances planetesimal-planetesimal 
collisions are even less important than at 2.5 AU (right panel of 
Figure |9j. The Earth-mass protoplanet has only a few parts in 
a thousand contribution from planetesimal collisions. 
Planetesimal growth at 25 AU 

In the outer regions of the protoplanetary disk we adopt a 
model where the surface density profile is proportional to the 
inverse of the orbital radius. This is in agreement with sur¬ 
face density profiles of observed protoplanetary disks (39, 40). 
Hence we set the total column density of particles at r = 25 
AU to Up = 0.54 g/cm 2 . We set the temperature at 25 AU 
according to full radiative transfer models of protoplanetary 
disks. This is important, since the optically thin assumption in 
the Minimum Mass Solar Nebula model of Hayashi (79) over¬ 
estimates the temperature in the outer disk. The growth rate of 
a planetesimal with mass M p by pebble accretion is given by 

-i 2 


Mp = 7r/i 


GM n 

(At 


(A v)p } 


P * 


(3) 


Here /b denotes the ratio of the actual accretion radius to the 
Bondi radius (24), Av is the sub-Keplerian speed difference 
and pp is the mid-plane particle density. The sub-Keplerian 
speed is given by 


Av = — 


1 H d\nP 

2 r 9 In r s 


(4) 


The speed difference scales with the sound speed as c 2 , since 
H/r contains another c s . Hence the accretion rate scales with 
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Figure 9: Growth of icy planetesimals at 10 AU. The growth rate by pebble accretion is as high as in the asteroid belt, since 
the lower column density of pebbles is counter-acted by the increased sedimentation in the more dilute gas. The annulus of 0.2 
AU width produces in the end an Earth-sized protoplanet and a single Moon-sized embryo. Pebble accretion overwhelmingly 
dominates the growth (right panel). The icy protoplanet that forms has only a few parts in a thousand mass contribution from 
collisions. Large Ceres-sized planetesimals have a contribution from collisions of less than 5%, while the 200 km planetesimal 
owes about 1/10 of its growth from 130 km to planetesimal collisions. Note how the Ceres-sized planetesimal (blue line in the 
right panel) got a head start for efficient accretion of pebbles by experiencing a significant collision after 700,000 years. 


the sound speed to the negative sixth power, and the sound 
speed must be set carefully to yield realistic results. We set 
the sound speed at 25 AU through the aspect ratio H/r = 0.05, 
based on radiative transfer models of protoplanetary accretion 
disks {41). This is much smaller than the standard value in the 
Minimum Mass Solar Nebula {H/r = 0.072 at r = 25 AU) 
and leads to much higher accretion rates in the outer nebula. 

We consider two models for planetesimal growth at 25 AU. 
First, a low density model where the internal density of the 
planetesimals is set to 0.5 g/cm 3 , similar to comets and some 
binary Kuiper belt objects (42). Second, a high density model 
that has instead an internal density of 2 g/cm 3 , more in agree¬ 
ment with the icy dwarf planets Pluto and Ceres. The low den¬ 
sity model experiences much lower pebble accretion rates be¬ 
cause of the reduced gravity of the planetesimals. We find that 
a turbulent stirring of a = 10 -4 , which we used for the aster¬ 
oid belt, results in growth times for pebble accretion of more 
than 10 Myr. Hence we set the turbulent stirring to a low value 
of a = 10 -6 in the low density model. The high density model 
has a = 10 -5 . We discuss in the next section how such low 
values could arise physically. 

We show the resulting size distributions and growth curve of 
the largest body at 25 AU in Figure [lO] The size distributions 
in the left panel of the figure both display a steep power law 
beyond 100 km in radius. Planetesimals up to 300 km in ra¬ 
dius accrete significant amounts of pebbles. Beyond 300 km a 


run-away growth sets in where the largest body detaches from 
the remaining population (right panel of Figure [TO}. This run¬ 
away growth is possible because of the weak gas friction at 25 
AU. Hence run-away pebble accretion is not slowed down by 
reaching the strong coupling branch but can instead continue 
to very large sizes (see Supplementary Material and Figures [S7] 
and[S3]>. 

Altogether it is clear that pebble accretion at 25 AU is much 
slower than at 2.5 AU, and actually requires very low turbu¬ 
lent activity for significant growth. But planetesimal growth 
in such wide orbits apparently also has the potential to result 
in run-away accretion up to ice giant sizes. The size distri¬ 
bution of planetesimals from 100 to 300 km in radius is very 
steep, in good qualitative agreement with what is observed in 
the present Kuiper belt (35). Our results show that 25 AU can 
be considered the very limit to where pebble accretion is sig¬ 
nificant, and then only for weak turbulence. Beyond this radius 
planetesimals maintain their birth sizes (and size distribution), 
and hence the very steep size distribution of large planetesimals 
in the cold component of the classical Kuiper belt may reflect 
the actual exponential cut off in the underlying birth size distri¬ 
bution. 

In the solar system, Neptune stopped its outwards migration 
at 30 AU. Therefore the primordial Kuiper belt can not have ex¬ 
tended much beyond that distance, as otherwise Neptune would 
have continued to migrate outwards by scattering planetesi- 
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Figure 10: Planetesimal growth at 25 AU. Two models are considered, a low density model where the internal density is set to 
p, = 0.5 g/cm 3 , similar to comets and binary Kuiper belt objects, and a high density model where the internal density is set to 
p, = 2 g/cm 3 , similar to the dwarf planet Ceres. The low density model has turbulent stirring a = 10” 6 while the high density 
model has a = 10 -5 . Both models display ordered growth up to 300 km radii, with a steep size distribution beyond 100 km 
sizes. This is followed by a run-away growth of a single, massive body. The right panel shows the size of the largest body as 
a function of time as well as the speed relative to a circular orbit. The run-away growth is facilitated by a steep decline in the 
eccentricity of the orbit, as the high pebble accretion rate damps the eccentricity. 


mals (38). However, the presence of the cold population, which 
appears detached from Neptune, indicates that the Kuiper belt 
could simply have transitioned to a much lower surface density 
around 30 AU. The planetesimals there did not grow beyond 
their birth sizes, resulting in a jump in the planetesimal column 
density between the regions inside 25 AU which experienced 
growth by pebble accretion and the regions outside which did 
not. 

Discussion 

Variation of turbulence strength 

An important and relatively unconstrained parameter in our 
chondrule accretion model for the formation of asteroids and 
planetary embryos is the strength of the turbulence in the pro¬ 
toplanetary disk, as the turbulent diffusion coefficient sets the 
scale-height of the particle layer and the mid-plane density of 
the chondrules (7). The asteroid formation region is located in 
a region of the protoplanetary disk where the degree of ionisa¬ 
tion is believed to be too low to sustain turbulence driven by the 
magnetorotational instability. The ionisation degree in the sur¬ 
face layers may nevertheless be high enough to drive accretion 
there. The turbulent motion in these layers will stir the oth¬ 
erwise laminar mid-plane. Effective turbulent viscosities be¬ 
tween a = 10 -5 and a = 10~ 4 have been measured in the 
mid-plane in computer simulations of dead zone stirring (43). 

Models of protoplanetary accretion disks that also include 


the effect of ambipolar diffusion in the surface layers find that 
the growth of the magnetorotational instability is quenched 
even there, by the lack of coupling between electrons and neu¬ 
trals (44). The magnetic field may enter a configuration where 
gas (and angular momentum) is expelled upwards, while gas 
connected to the magnetic field lines closer to the mid-plane is 
accreted towards the star. This disk wind accretion leaves the 
mid-plane completely laminar. The sedimentation of dust will 
nevertheless cause a mild stirring of the mid-plane (45). The 
stirring caused by streaming and Kelvin-Helmholtz instabili¬ 
ties in the dense mid-plane layer of chondrules will be much 
milder than the stirring from the active layers, with a typical 
value of a = 1(T 6 for chondrule-sized particles. 

The nominal value for the turbulent diffusion coefficient in 
our model is a = 10 -4 , corresponding to conditions in a dead 
zone stirred by active surface layers. To quantify the effect 
of the turbulent stirring of chondrules on our results we have 
run simulations with a lower (a = 2 x 1CV 6 ) and a higher 
(a = 10 3 ) value of the turbulent diffusion coefficient. The 
size of the largest planetesimal in the asteroid belt is shown as 
a function of time in Figure [XT] Clearly, a lower value of the 
turbulent diffusion coefficient results in a much higher growth 
rate of the planetesimals, as chondrules are allowed to sediment 
much more strongly to the mid-plane. For a = 10 -3 the for¬ 
mation of the first embryos is delayed to after 5 Myr, which is 
longer than the typical life-times of protoplanetary disks around 
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Figure 11: The maximum planetesimal radius in the asteroid 
belt versus time, for three different values of the turbulent vis¬ 
cosity a. Here a = 2 x 1(T 6 represents the strength of turbu¬ 
lence caused by streaming instabilities and Kelvin-Helmholtz 
instabilities in a sedimented mid-plane layer of chondrules, 
a = 10' 1 represents the turbulence strength in a dead zone 
stirred by active surface layers and a = 10~ 3 the turbulence 
strength caused directly by the magnetorotational instability. 
Turbulent stirring of chondrules sets the scale-height and mid¬ 
plane density of the chondrule layer and hence dictates the 
planetesimal growth rate. The formation time of the first em¬ 
bryo depends strongly on the degree of stirring. 

young stars. 

Layered accretion in the asteroid belt 

Our model provides a direct connection between chondrules 
and asteroid growth. While early accreted asteroids and their 
chondrules would melt from the energy released by the decay 
of the short-lived 26 A1 radionuclide, layers accreted later than 
2 Myr may preserve their pristine, undifferentiated nature (46). 
Indeed, the existence of differentiated asteroids overlaid by 
chondritic crusts has been proposed to explain the systematic 
magnetisations of chondrules in the Allende meteorite (47,48), 
imprinted from a geodynamo operating in the liquid central re¬ 
gion of the parent body. The asteroid Lutetia has been sug¬ 
gested to be partially differentiated, due to its high mean den¬ 
sity and primitive surface (49). 

The characteristic maximum size of planetesimals formed 
through streaming instabilities is 50-150 km in radius for col¬ 
umn densities comparable to the Minimum Mass Solar Nebula 
(Figure [3]> . Planetesimals much smaller than 100 km in radius 
are relatively inefficient at accreting chondrules and hence ac¬ 
crete only a thin veneer of chondrules on top of the primor- 
dially formed planetesimal (Figure [5] and Figure |5J. Planetesi¬ 


mals forming before 1-2 Myr after CAIs will likely differenti¬ 
ate, as was the case for the parent bodies of the differentiated 
meteorites (50, 51). Hence we predict that asteroids currently 
residing in the knee at the size distribution at 60 km radii rep¬ 
resent primordial, differentiated bodies overlain by a veneer of 
chondrules of thickness up to a few dozen kilometers. 

Asteroids larger than 100 km in radius, on the other hand, 
grow mainly by accretion of chondrules. The initial planetes¬ 
imal seeds will be covered with extended layers of narrowly 
size-sorted chondrules. The Allende meteorite, whose system¬ 
atic magnetisation makes the parent body a prime candidate for 
layered accretion, has been proposed to originate from the Eos 
family in the asteroid belt (52). This family is known to display 
significant spectral variations between its members (52), indi¬ 
cating partial differentiation and internal heterogeneity of the 
110-km-radius parent body. Another example of an asteroid 
family with large internal variation in the spectral properties is 
the Eunomia family (48), believed to originate from a parent 
body of 150 km in radius. Intermediate-sized asteroids of 100- 
200 km in radius appear to be the best candidates for produc¬ 
ing such families with large internal variation. Bodies of this 
size range have comparable mass fractions in the differentiated 
planetesimal seed and in the accreted chondrule layers. 

Many other asteroid families display internal albedo distribu¬ 
tions that are much narrower than the spread of albedos across 
families (55), in stark contrast to the heterogenous Eos and Eu¬ 
nomia families discussed above. The inference of the internal 
structure of the parent bodies of such apparently homogeneous 
families is nevertheless complicated by the identification of in¬ 
terlopers in (and exclusion from) the asteroid families. A fam¬ 
ily member whose albedo is different from the family mean 
could in fact be mistaken for an interloper. However, the iden¬ 
tification of asteroid families that appear to originate from inter¬ 
nally homogeneous asteroids is not in conflict with our model. 
Our results show that many asteroids in the 50-100 km radius 
range will have experienced little chondrule accretion and rep¬ 
resent primordial asteroid seeds. Additionally, very large as¬ 
teroids of radii larger than 300^400 km must have undergone 
significant accretion of chondrules between 2 and 3 Myr after 
the formation of the seed planetesimal. Families formed mainly 
from the chondrule layers of such large asteroids would also ap¬ 
pear relatively homogeneous. Fragments as large as 100 km in 
radius could consist of purely chondritic material, and families 
produced from those fragments would in turn appear entirely 
homogeneous. 

Implications for terrestrial planet formation 

The largest planetesimals in our simulations of the asteroid 
belt and the terrestrial planet formation region grow to sizes of 
several thousand kilometers, forming the embryos whose mu¬ 
tual collisions drive the subsequent terrestrial planet formation 
stage. With only 11% of Earth’s mass. Mars may be one of 
these remaining embryos. Mars is inferred to have accreted 
2 ± 1 Myr after the first condensations in the solar system (54), 
in good agreement with the formation time-scale of planetary 
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embryos by chondrule accretion in our simulations. This im¬ 
plies that chondrule accretion is a driving mechanism ensuring 
the rapid growth of planetary embryos in the terrestrial planet 
formation region. Thus the production of chondrules may be a 
key process promoting the assembly of rocky planets. 

Materials and Methods 


Planetesimal formation simulations 

The planetesimal formation simulations were per¬ 
formed with the Pencil Code (which can be freely down¬ 
loaded, including modifications done for this work, at 
http://code.google.com/pencil-code). A new sink particle 
algorithm developed for this project is presented in the 
Supplementary Material. 

Our simulations of planetesimal formation through stream¬ 
ing instabilities are performed in the local shearing box frame. 
The simulation box orbits with the local Keplerian frequency 
17 at an arbitrary radial distance r from the central star. We fix 
the box size to a cube of side lengths L = 0.2 H, where H is 
the gas scale-height. The box size is chosen in order to cap¬ 
ture the relevant wavelengths of the streaming instability (55). 
The gas is initialised with a Gaussian density profile around 
the mid-plane, in reaction to the component of stellar gravity 
directed towards the mid-plane. The temperature is assumed 
to be a constant set through the sound speed c s , which is kept 
fixed during the simulation. Particles are given positions ac¬ 
cording to a Gaussian density distribution, with a scale-height 
of 1% of the gas scale-height. We assume that particles in the 
entire column have sedimented into the simulation domain, giv¬ 
ing a mean solids-to-gas ratio of (p p ) / p 0 ss 0.25 for the chosen 
metallicity Z = 0.01. The particles have uniform sizes given 
by the Stokes number St = 17/f = 0.3, where Tf is the friction 
time of the particles, corresponding to approximately 25-cm- 
sized particles at 2.5 AU. It has previously been shown that the 
consideration of a range of particles sizes gives similar results 
to the monodisperse case (12). 

The strength of the particle self-gravity in a scale-free local 
box simulation is set by the non-dimensional parameter 
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Here G is the gravity constant, p (l is the gas density in the 
mid-plane and f2 is the Keplerian frequency at the chosen ra¬ 
dial distance from the star. The scaling with semi-major axis 
is based on the properties of the Minimum Mass Solar Neb¬ 
ula (79). The parameter r multiplied by the non-dimensional 
particle density p p /po enters the Poisson equation for self¬ 
gravity, V 2( P = AnGp p , when lengths are measured in units of 
the scale-height and times in units of the inverse Keplerian fre¬ 
quency f2~ - 1 . The parameter also converts a given mass density 
of particles in a grid cell to an actual mass, since the choice of r 
defines the mid-plane gas density po. Thus both the dynamics 
and resulting planetesimal masses are strongly influenced by 
the value of /’ (this is also evident from the results presented in 
Figure |5J. 


Chondrule accretion simulations 

We have developed a statistical code which evolves the 
masses and orbits of planetesimals as they accrete chondrules 
and collide mutually. Details of the code algorithms, test cal¬ 
culations as well as comparisons to other published results are 
presented in the Supplementary Materials and Methods. 
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Supplementary Materials 

A. Johansen, M.-M. Mac Low, P. Lacerda, & M. Bizzarre 


1 Materials and Methods 

1.1 Sink particles in the Pencil Code 

In order to track the masses of the planetesimals as they form and grow, we have developed a new sink particle module for the 
Pencil Code. Above a threshold particle density value p p = pp Smk) all particles in the cell are merged to a single sink particle. 
This sink particle does not feel friction with the gas, due to the large size of the planetesimal. Any particle coming within the 
distance 5x/2 of a sink particle is destroyed (Sx is the size of the resolution element) and its mass and momentum added to 
the sink particle. In Figure [sT] we compare the planetesimal masses arising from simulations at 256 3 with three different values 
for the sink particle creation threshold. The figure shows that the results are relatively unaffected by the choice of sink particle 
creation thresholds above 5 times the Roche density. Sink particles allow us to easily keep track of the temporal evolution and 
saturation of the planetesimal masses, and the merging of several million superparticles into a few sink particles prevents the 
simulation from slowing down when superparticles cluster on a subset of the available processors (56). 

1.2 Chondrule accretion simulations 

While classical simulations of planetesimal growth focus on the growth by gravitationally focused planetesimal-planetesimal 
collisions (20, 21, 57), the accretion of chondrule-sized objects coupled to the gas via friction can potentially be much more 
relevant if planetesimals are born and grow in an ocean of chondrules (23,24,58). “Chondrule accretion” is an aspect of the more 
general term “pebble accretion”, but chondrules are much smaller than the typically cm-sized pebbles used in previous studies of 
pebble accretion. 


1.3 Pebble accretion 

Pebble accretion has two dominant regimes. Planetesimals below a transition radius of around 1000 km (at 2.5 AU) accrete 
pebbles embedded in the sub-Keplerian gas. In the Minimum Mass Solar Nebula (19) the relative speed between a planetesimal 
on a circular orbit and the gas is Aw ~ 50 m/s. These pebbles are accreted by the planetesimal from the Bondi radius of the 
planetesimal, 

Rb = GM/(Av) 2 , (6) 

provided that the friction time of the pebble matches the time to cross the Bondi radius. Larger pebbles are simply scattered by 
the planetesimal, while smaller pebbles can not be pulled away quickly enough from the almost straight streamlines of the gas. 
The finite size of the planetesimal leads to additional complexities, namely gravitational focusing of loosely coupled pebbles; 
and the Stokes flow of the gas around the planetesimal, which affects the trajectories of strongly coupled particles with impact 
parameters similar to or smaller than the radius of the planetesimal. 

In order to model all these effects correctly we have solved for the trajectories of a large number of pebble-planetesimal 
combinations, for different values of the friction time tf and the planetesimal radius R. The planetesimal is kept stationary at 
(x, y ) = (0,0), while the pebble moves under the influence of planetesimal gravity and gas drag. Pebbles enter the simulation 
domain at a large, positive value of x and with impact parameter b = y 0 . The gas velocity field is kept fixed and follows the 
Stokes solution (59) 


u r 


Ug 


—Aw cos (6) 
+Aw sin(0) 


3 R 1 R 3 \ 


2 r 2 r 3 J 

(7) 

3 R 1 R 3 \ 


4 r 4 r 3 / 

(8) 


Here r and 0 are the polar coordinates, with r denoting the distance from the origin and 6 the angle between the rr-axis and 
the y-axis, so that the gas flow at large r is uniform along the ^-direction u = —Avx. Collisions with the planetesimal are 
treated as instantaneous, conserving the total momentum as well as the speed component parallel to the surface and a fraction c 
of the speed component perpendicular to the surface, with c denoting the coefficient of restitution. The pebble is assumed to be 
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Figure SI: The non-dimensionless planetesimal masses // = p p (Sx) 3 (here p p is the particle density represented by the planetesi- 
mal in a grid cell and Sx is the size of the cell) and corresponding contracted radii as a function of time after self-gravity is turned 
on. Simulations at 256 3 grid cells and three different values of the sink particle creation threshold are compared. Creation of a 
planetesimal is marked with an empty circle, while the destruction is marked with a filled circle. The sizes of the largest plan- 
etesimals are relatively unaffected by the choice of the threshold. There are more mergers of small planetesimals at a low sink 
particle creation threshold, but the planetesimals resulting from those collisions correspond well to the massive planetesimals 
forming at higher creation thresholds. 


accreted after colliding with the planetesimal 5 times or after orbiting the planetesimal 5 times. The latter criterion is convenient 
for planetesimals with large Bondi radii where the pebble enters a slowly decaying orbit around the planetesimal. 

We show the accretion radius R acc (the impact parameter required for accretion) for (a) perfect sticking and (b) a coefficient 
of restitution of c = 0.5 in Figure [S2| The accretion radius depends on both the friction time (normalised by the Bondi time 
t| S = /i'| S /Au in Figure [S2] > and the planetesimal size (normalised by the Bondi radius). We discuss the general trends in the 
accretion radius curve here: 


Large Bondi radius. The gas flow around the planetesimal and the coefficient of restitution play no role in determining the 
accretion radius when the planetesimal is much smaller than the Bondi radius. This is the case for large planetesimals and/or 
low relative pebble-planetesimal speeds. The three sections on the growth curve in Fig. S2 are: (i) Strong coupling when the 


pebble couples to the gas on a much shorter time-scale than the Bondi time-scale. Here the accretion radius is proportional to 
Rb s/tf/tB- (ii) Weak coupling when the pebble couples to the gas on much longer time-scales than the Bondi time-scale. Here 
the accretion radius drops rapidly with increasing friction time, with no simple analytical expression, (iii) Gravitational focussing 
for the largest particles when the weak coupling branch drops below the gravitational focussing radius. 


Small Bondi radius. A small planetesimal has Bondi radius much smaller than the planetesimal radius. This limit has two 
distinct branches: (i) Sedimentation overshoot for strongly coupled pebbles. These pebbles follow the Stokes flow around the 
planetesimal, but manage to sediment on to the surface for very low impact parameters where the gas flow comes very close to 
the planetesimal surface, (ii) Geometric accretion for weakly coupled pebbles that do not react to the modified gas flow close to 
the planetesimals. The latter branch in turn is strongly affected by the coefficient of restitution. Perfect sticking leads to accretion 
onto the entire projected surface, while consideration of the collision outcome requires very low impact parameters to remain 
bound after the collision. 


Below we discuss some important aspects of the sedimentation overshoot and geometric branches. 

Sedimentation overshoot. Small particles are carried with the gas around the planetesimals. The frictional acceleration between 
gas and particles is nevertheless not instantaneous but operates on the time-scale t{. This can lead to accretion of particles 
that sediment onto the planetesimal surface during the transport around the planetesimal. The gas streamlines at a low impact 
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Figure S2: The accretion radius i? acc (i.e. the impact parameter required for accretion, given in units of the Bondi radius 
versus the particle friction time if (in units of the Bondi time-scale t B ), for planetesimal radii R between 10 -4 f?B and 10 4 R B . 
The left panel shows the results when assuming perfect sticking between pebble and planetesimal. Planetesimals with R < Rb 
have three rather distinct branches: strong coupling for short friction times, weak coupling for longer friction times and finally 
gravitational focusing for the longest friction times. Planetesimals with R > Rb radius have two branches: sedimentation 
overshoot for short friction times (which follows the strong coupling scaling but is physically distinct) and geometric accretion 
for friction times longer than the time to pass the planetesimal. The consideration of the outcome of the pebble-planetesimal 
collision, with a coefficient of restitution c = 0.5, is shown the right panel. This strongly limits geometric accretion onto small 
planetesimals. 


parameter of yo = b <C R move around the planetesimal at the distance Sr ~ b, at the approximate speed 


u e ~ Au— . 

Ignoring factors of order unity, the time to move around the planetesimal is then 

R R 2 
e us bAv 

During this time the pebble sediments towards the planetesimal at the terminal velocity 

GM 

The criterion for reaching the surface is vt tq ~ b. This in turn yields 

,2 GM 
b 2 ~ tf 


Av 


Dividing by the squared Bondi radius finally yields the criterion for accretion 

b 2 tf 

Rb t B ' 


(9) 


( 10 ) 


( 11 ) 


( 12 ) 


(13) 


This expression is equivalent to the strong coupling limit of pebble accretion (24), although the physics of the accretion via 
sedimentation overshoot is completely different. 
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Figure S3: The evolution of eccentricity e and inclination i of 1000 planetesimals with mass M = 10 24 g located at 1 AU with a 
surface density of 10 g/cm 2 , for three values of the time-step. Here Cdt denotes the time-step relative to the collision time-scale. 
Stirring and dynamical friction are implemented using the scheme described in Ohtsuki et al. (2002). The results are qualitatively 
similar to Figure 4 of that paper, although the e and i evolve a bit slower in our simulation. We attribute this to a difference in the 
set up of the simulations: while we initially set e = 10 -4 and i = 5 x 10 1 for all planetesimals, Ohtsuki et al. (2002) give their 
planetesimals a Rayleigh distribution around those values. We use cat = 0.1 in the actual planetesimal growth simulations. 


Geometric accretion. The requirement that the pebble remains bound after the collision can be found geometrically as 


Rare. 


{v e /Av) 2 — C 2 
1 — C 2 


(14) 


Pebbles entering with an impact parameter less than R :iCC are thus accreted, while pebbles with a larger impact parameter collide 
once with the planetesimal and then fly off to infinity. The accretion radius is only positive for c < v e /Av = \/‘2 Ry>J R. Hence 
planetesimals with Bondi radii smaller than their size can not accrete pebbles which couple to the gas on time-scales longer than 
the time-scale to pass the planetesimal (unless c is extremely low). 


1.4 Eccentricity and inclination 

The growing asteroids experience mutual gravitational encounters that excite the eccentricities and inclinations of the population. 
This affects the pebble accretion rate since the relative speed between a planetesimal and the sub-Keplerian pebble flow changes 
with position in the orbit for non-zero eccentricities and inclinations. 

The temporal evolution of the eccentricity e and the inclination i is calculated using the analytical approximations of Ohtsuki 
et al. (52), constructed to match the results of TV-body simulations for both low and high eccentricities. The analytical evolution 
equations for e and i give the excitation as well as dynamical friction of a binned population of planetesimals. We show in Figure 
|S3|the result of a test problem defined in Ohtsuki et al. Here we consider 1000 equal-mass planetesimals with initial eccentricity 
of 10 -4 and inclination of 5 x 10 -5 . The evolution of e and i follows the general curve in Figure 4 of Ohtsuki et al., but the 
evolution is a bit slower. We attribute this to a difference in the tests. In Ohtsuki et al. the planetesimals are given a Rayleigh 
distribution of e and i initially, while we assign a constant value equal to the mean of that Rayleigh distribution. A similar test 
problem with 800 planetesimals is shown in Figure[S4] with good agreement with Figure 2 of Stewart & Ida (60). 

In principle we could use all the planetesimal particles in the code as natural bins and evolve each planetesimal with the 
dynamical excitation contribution from all other planetesimals. However, this would take unfeasibly long time to compute when 
there are millions of planetesimal particles. Instead we bin the planetesimals by their radius and consider the dynamical evolution 
of the smallest and largest planetesimal in each bin. Gravitational stirring is only considered between these anchor planetesimals. 
The resulting values of de 2 /df and dt 2 /df are then multiplied by the total number of planetesimals in the stirring bin, divided 
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Figure S4: The evolution of eccentricity e and inclination % of 800 planetesimals with mass M = 10 24 g located at 1 AU with a 
surface density of 10 g/cm 2 . The eccentricity and the inclination are normalised by their Hill values [2M p /(3M*)] 1//3 . Results 
are qualitatively similar to Figure 2 of Stewart & Ida (2000). As in Figure [S3] a time-step parameter of c: dt = 0.1 provides a 
good compromise between precision and speed; hence we use c<jt = 0.1 for planetesimal stirring in the planetesimal growth 
simulations. 


by the number of planetesimals that were actually considered (the two anchors). The evolution of e and i is finally interpolated 
from the anchor planetesimals to all the other planetesimals at the end of the time-step. 

The eccentricity and inclination of the orbit can now be used to calculate the pebble accretion rate over a full orbit of the 
planetesimal. For small eccentricities and inclinations the orbit can be considered in the local coordinate frame corotating with 
the Keplerian flow at a given distance ro from the star. The coordinate system is oriented such that the .x‘-axis points radially 
away from the star, the y-axis along the orbital direction of the gas and the ::-axis perpendicular to the plane of the disk, along 
the rotation vector of the disk. The epicyclic motion of a planetesimal with eccentricity e in such a frame is 


V x (t) 

= W e C0s(f2f) , 

(15) 

V y (t) 

= -^w e sin(f?f), 

(16) 

v z (t) 

= WiCOS (fit). 

(17) 


Here w e = ev K () and W; = Wk,o are the eccentricity and inclination speeds of the orbit, respectively, measured relative to the 
Keplerian speed at the centre of the frame, and ( v x ,v y ) is the velocity is measured relative to the local Keplerian motion at the 
instantaneous position of the planetesimal. Relative to wk,o the azimuthal velocity is 

v y (t) = — 2w e sin(f?f). (18) 

The velocity of the incoming chondrules is v' y = — Aw relative to the local Keplerian velocity. Hence the relative speed between 
planetesimal and chondrules is 


w re i = \J [w e cos(fit)] 2 + [—(l/2)w e sin( fit) + Aw] 2 + [w; cos(l?t)] 2 . (19) 

The relative speed is lowest at aphelion (fit = tt/2) and highest at perihelion (fit = 3 tt/2). Since the Bondi radius scales as 
Rb oc l/w 2 el , the eccentric motion strongly affects the accretion rate. On the strong coupling branch of pebble accretion one can 
nevertheless show that the eccentric orbit will not affect the mass accretion rate, because in fact the accretion rate is independent 

of w re i, 

t £ Q2 M 2 

M sc oc —Rb Wrei OC 3 -4 —w re i oc GMtf . (20) 

t B GM/w 3 el w r 4 eI 
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Figure S5: Phase points covering planetesimal orbits with inclinations i = 0, i = 0.0016, i = 0.0032 and i = 0.0064. We 
construct phase arrays to cover both the full planetesimal orbit and the scale-height of the chondrule layer (here H p = 0.02 H). 
The height in the planetesimal orbit is normalised by the maximum height for the given inclination. The dotted line indicates the 
scale-height of the chondrule mid-plane layer, clearly resolved for all choices of the inclination. 


At aphelion the relative speed between chondrule and planetesimal can even approach zero for (l/2)w e ~ Aw. At this point the 
Bondi radius is no longer the relevant accretion radius as the Hill speed of the planetesimal is higher than the relative chondrule- 
planetesimal speed, and the planetesimal thus enters short periods of Hill accretion at aphelion (our implementation of Hill 
accretion is explained in section [L6| ). 

We take into account the eccentric and inclined orbit of the planetesimal by sampling the accretion rate at a number of phases 
in the orbit. The choice of phase points is a balance between the need to sample both the scale height H p of the chondrule layer, 
where chondrule densities are highest, as well as the planetesimal orbit over which the relative planetesimal-chondrule speed, 
and hence the accretion radius, varies. We set the first phase point at $ = = 0. The distance to the next phase point is taken 

as the minimum of the two functions 


A$i 

A $ 2 


exp(z/H c )H c 27r 

A>rb A| 

27T 


n 2 ■ 


( 21 ) 

( 22 ) 


Here z OI b is the maximal height of the planetesimal relative to the mid-plane. This gives phases that cover the full planetesimal 
orbit, with additional resolution elements added within the chondrule scale height. This way we avoid having to use a very large 
number of points for highly inclined orbits. We show examples of phase arrays in Figure[S5]for N\ = 48 and N 2 = 12. 

In Figure [S6] we show the mass accretion rate of planetesimals at 2.5 AU as a function of the eccentricity of the orbit. We use 
the approximation that i = e/2 and measure the accretion rate relative to that at e = i = 0. Increasing the eccentricity actually 
increases the geometric accretion rate on small planetesimals (below 50 km in radius), since the relative planetesimal-chondrule 
speed increases and hence the mass flux increases proportionally to the speed. Larger planetesimals are affected negatively by 
eccentricity. These planetesimals accrete small chondrules at very high rate on circular orbits. The accretion rate at the aphelion 
of an eccentric orbit does not benefit strongly from the decreased relative planetesimal-chondrule speed at that phase, since the 
larger Bondi radius requires very large chondrules, which are not present in the disk. Instead the planetesimal enters the strong 
coupling branch at aphelion, while the accretion rate is strongly reduced at perihelion. The overall result is a reduction in the 
accretion rate. Even larger planetesimals are affected less and less by the increase in eccentricity. These large planetesimals 
already accrete on the strong coupling branch, which is relatively unaffected by the relative speed and hence by the eccentricity. 

Figure[S7]illustrates the mass accretion rates of planetesimals at 25 AU for different values of their eccentricity. The normalised 
accretion rates in Figure[S7]are flat or increasing beyond 200 km in radius, indicating a very steep increase in the accretion rate 
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Figure S6: The effect of eccentricity and inclination on the chondrule accretion rate, at an orbital distance of 2.5 AU. The x-axis 
shows the planetesimal sizes and the y -axis the pebble accretion rate normalised to the accretion rate from the full Bondi radius. 
Dotted lines mark 50 km, 100 km and 200 km. We use here the approximation i = el 2. Increasing the orbital eccentricity 
can actually increase the accretion rate on bodies smaller than 100 km in radius, as the orbital speed in aphelion matches the 
sub-Keplerian speed of the chondrules which leads to very high accretion rates. Larger planetesimals are affected negatively 
beyond the threshold (e/2)uK = At; as the large particles necessary for efficient accretion at aphelion are not present in the disk. 
The top plot shows results for the column density of the Minimum Mass Solar Nebula at 2.5 AU, the middle plot 0.3 times the 
MMSN and the bottom plot 0.1 times the MMSN. 


with size. This is in contrast to the situation in the asteroid belt where the run-away accretion of chondrules is stopped by the 
lack of chondrules larger than mm in size (Figure [S6]). 
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Figure S7: The accretion rate of planetesimals at 25 AU, as a function of the planetesimal size. The accretion rate is normalised 
by the Bondi rate. The top panel shows the accretion rate for pebbles up to 1 mm in radius, while the bottom panels shows the 
accretion rate for pebbles up to 1 cm in radius. The accretion rate beyond 200 km in radius is flat or increasing here in the outer 
regions of the protoplanetary disk. This is in contrast to the situation in the asteroid belt where the lack of chondrules of larger 
than mm size slows down the run-away accretion of large planetesimals (Figure [S6|. The presence of mm-cm pebbles in the outer 
protoplanetary disk can drive a run-away accretion where the most massive planetesimals detach from the rest of the population 
(see Figure p~0]>. 

1.5 Stratification integral 

Calculation of the mass accretion rate requires knowledge of the accretion radius as well as the chondrule density averaged over 
the accretion radius. The stratification integral S is defined as the mean chondrule density normalised by the chondrule density in 
the mid-plane. The stratification integral for a planetesimal with accretion radius R acc located at the height z 0 over the mid-plane 
is 

| /*^0 - l - -^acc 

S =—^T ex p[—2 2 /(2-ffp)]2i/ Rl cc -(z- z 0 ) 2 dz . (23) 

7r ' K acc Jzo — Racc 

This expression is obtained by summing over lines of constant z and hence constant chondrule density. There is no simple 
analytical solution to equation ( |23j ). Tabularisation of the numerical solution requires interpolation in both f? acc , zq and the 
angle of incidence of the chondrule flow 0 (see below). We therefore an approximation that allows the integral to be calculated 
analytically. This square approximation integrates the chondrule density over a square instead of a circle. This yields the solvable 
integral 

i /*zo+rt acc 

Square = , / exp[-z 2 /(2H^)}2R acc dz . (24) 

^-^aCC J Zq — Ra, cc 
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Figure S8: The stratification integral (i.e., the mean particle density over the accretion cross section) as a function of the accretion 
radius of the planetesimal. Numerical integration of the Gaussian stratification over a circle (blue lines) is expensive, so we use 
instead a square approximation to the integral (red lines). Full lines indicate a planetesimal in the mid-plane, dashed lines a 
planetesimal at two times the chondrule scale-height above the mid-plane, and dash-dotted lines a planetesimal at ten times the 
chondrule scale-height. The square approximation is very precise both when the accretion radius is much lower or much higher 
than the chondrule scale-height. 


The analytical solution is 


, _ H p v/tt 

' square V2R acc 2 


erf 


( Z0 + Race | 

V V2H P ) 


— erf 


| ZO Race \ 

V V2H p ) 


(25) 


The square approximation tends to underestimate the chondrule density because of the inclusion of low-density corners in the 
square. Therefore we decrease the size of the square by replacing f? acc in equation ( |25[ ) with R' acc = /f? acc ■ We have found 
/ = 0.79 to give a much better fit to the numerical integral than / = 1. In Figure [S8| we show how the square approximation 
compares to the full numerical integration of the accretion radius over the Gaussian stratification. 

Planetesimals at very high inclinations encounter the chondrule flow at an angle 8 which is / 0°, measured relative to the 
vertical. The angle of incidence is found through 

/ (re 1)\ 

8 = acos ( —— I (26) 

y I'rel J 

Here ri re i = — v c \ is the vectorial relative speed between the planetesimal (with velocity v p ) and the chondrules 

(with velocity v c — —A vy). The components of v p are given in equations ([T5])-<[T7]). The consideration of the angle of incidence 
changes the stratification integral to 


pZ o+-Racc cos 0 


S = 


nRlcc jZo-Rz.cc cose 


exp[-z 2 /(2H 2 )]2\ Rl cc - 


z — Zo \ “ dz 
cos 9 ) cost? ’ 


(27) 


The square approximation changes to 


5 S , 


H p y/n 1 
V2Ro.cc 2 cos( 


e Zo + Race COS 8 \ ( Z 0 - i? acc COS 8 
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The results presented in this paper are obtained with the square approximation including the correction for the encounter angle. 
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1.6 Transition to Hill accretion 

Bondi accretion is only valid for planetesimals below a transition mass of approximately 0.1% of an Earth mass (24). Beyond 
this the Hill radius of the growing planetesimal, /in. is so large that the relative speed between chondrule and planetesimal is set 
by the Hill speed Uh = fi/in rather than the sub-Keplerian speed Av. 

We obtain a smooth transition from Bondi accretion to Hill accretion by solving for the combined speed v(R acc ) = Av + 
f2R acc . The accretion radius R acc depends on the approach speed, so we solve the equation iteratively from the starting point 
v = Av. This allows us to also apply the tabulated accretion radius for Bondi accretion (Figure |S2[ ) to Hill accretion, by 
modifying the approach speed to take into account the Keplerian shear. We have checked that the resulting accretion rates on the 
Hill branch are very similar to the results from hydrodynamical simulations (24). 


1.7 Eccentricity damping by gas drag and chondrule accretion 


The planetesimals in our model are so large that 
relative speed squared and with friction time 


the quadratic drag force regime applies, with drag force proportional to the 


tf 


6 Rp, 
6vp g 


(29) 


The drag force can be applied directly to the eccentricity and inclination, following Wetherill & Stewart (1989) (61) and Mor- 
bidelli et al. (2009) (27), 


de 2 
d t 
di 2 
d t 


16 v Q.5irp g v 2 R 2 
l 42 M(l + 0.8/3 2 ) ’ 

o oo 2 v 0.5np g v 2 R 2 
' P 4- 2M(1 + 0.8/3 2 ) ' 


(30) 

(31) 


Here we use the notation of Morbidelli et al. (2009) who define (3 = i/e, v = fK\/(5/8)e 2 + (l/2)i 2 and v g = \Jv(v + At;). 
This formulation damps the eccentricity and inclination to zero on the relevant friction time-scale, with the term v + Av in v g 
corresponding to an orbitally averaged relative speed which would enter equation ( [29] ). We have checked the analytical evolution 
of e and i against an 7V-body integration of an eccentric orbit circularised by gas drag and found excellent agreement. 

Chondrule accretion and scattering are two additional damping mechanisms. The friction time-scale for chondrule accretion 
with accretion radius R acc is 


7f,cho — 


M 


4 R 2 


Pg 


nR 2 acc P v 5v 18 7? 2 CC p c 




gas • 


(32) 


We use similar expressions to equations (30i-(31 1 for the damping by chondrule accretion, but with R 1 replaced by R: 


2 

acc’ 


p& 


replaced by p p , and multiplied by 18/4 to take into account that the pre-factor in the friction time of chondrule accretion (4/3) 
is much lower than for gas drag (6). Damping by bouncing collisions and chondrule scattering is taken into account by using the 
maximum of the Bondi radius and the physical radius as the gravitational interaction radius when tf > ig and R acc < R- 


1.8 Planetesimal growth simulations 

The dynamical equations describing the temporal evolution of the masses, eccentricities and inclinations of a large number of 
planetesimals are solved in a numerical code named “Pebble Accretion Onto Planetesimals and Planets” (PAOPAP), which we 
have developed for the purpose of demonstrating the importance of chondrule accretion for planetesimal growth. Information 
about the simulations can be also found in the main paper; here we give some additional details. 

The planetesimals are treated as individual particles in the code, each represented by a mass, an eccentricity and an inclination. 
Planetesimal masses are evolved via the accretion of chondmles and mutual planetesimal collisions (see Section |T~9| ), while the 
eccentricities and inclinations are changed by both viscous stirring and dynamical friction from the other planetesimals, as well 
as damping by gas drag and pebble accretion and scattering. We ignore mass erosion of small (< 50 km) planetesimals by large 
chondrules impacting and bouncing at super-escape speeds. 

The code divides the planetesimals into discrete size bins and calculates M, e and % for the smallest and largest planetesimal 
in each bin. The temporal evolution of all other planetesimals is found by interpolation from the two anchors in each bin. This 
approach allows us to simulate a high number of planetesimals at a relatively low computational cost and to resolve the run-away 
growth of a small number of large objects. We use 200 logarithmically spaced bins spanning from 10 km to 10,000 km in radius. 
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The chondrule density has a Gaussian stratification profile, with the scale-height H p in each bin set according to the diffusion- 
sedimentation equilibrium expression 

H e 


fit{ + a 


(33) 


Here H g is the gas scale-height, a is the turbulent viscosity (which we assume takes the same value as the turbulent diffusion 
coefficient), 17 is the Keplerian frequency at the orbital distance to the asteroid belt and if is the friction time. The latter is given 
in the Epstein regime by 


tf 


ap. 


(34) 


CsPg 


with a denoting the particle radius, p, the material density, c s the sound speed in the gas (which depends only on the temperature 
and the mean molecular weight), and p g the gas density. We choose a nominal turbulent diffusion coefficient of a = HP 1 
(we discuss this choice in the main paper). This yields a scale-height of mm-sized particles of approximately 30% of the gas 
scale-height, with smaller particles having scale-heights increasing as the square root of their size. 

We create chondrules continuously on a time-scale of 1.5 Myr to satisfy the cosmochemical evidence that chondrules in 
individual chondrites have ages from 0 to 3 Myr (5). We have performed test simulations where all the chondrules are present 
from the beginning; this accelerates chondrule accretion but does not affect the overall picture that asteroids, Kuiper belt objects 
and planetary embryos grow by chondrule accretion. The chondrules are divided into 30 logarithmically spaced bins, with the 
number density n(a) distributed according to the size distribution dn(a)/da oc a ~ 3 5 . 

The planetesimals have initial sizes from 10 to 150 km, distributed in a shallow power law with the differential number 
following dN/dR oc R~ 2 8 , but truncated with an exponential term exp[—(7?/7? exp ) 4 ] at the planetesimal radius f? exp =100 km. 
The smallest planetesimals experience very little growth and act mainly to provide dynamical friction, damping the inclinations 
and eccentricities of the larger asteroids and planetary embryos. 

The pebble accretion rates of the planetesimals are found via interpolation in a look-up table that gives the accretion radius for a 
grid of values of the planetesimal size R/Rb (normalised by the Bondi radius) and the chondrule friction time if/fg (normalised 
by the Bondi time). The look-up table is based on a high number of integrations of the dynamics of single chondrules passing by 
a planetesimal with the sub-Keplerian flow. 

The temporal integration of the dynamical equations is done via a simple Eulerian scheme, with the time-step determined to 
make M, e and i change by a maximum of 10%. We experimented with 30% and found significant changes to the results, while 
3% gives very similar results. 


1.9 Planetesimal collisions 


Planetesimal collisions are included in the code via a Monte Carlo method. The planetesimal particles are first sorted in discrete 
size bins. For each bin we calculate the average inclination and eccentricity of the planetesimals. We then calculate the collision 
rate matrix for all combinations of bins, i\j. The collision rates are calculated following the scheme described in the online 
supplement of Morbidelli et al. (2009) (21). We refrain from describing the scheme in detail here since it is already well 
described in Section 1.1 of the online supplement of Morbidelli et al. (2009) (27). 

From the time-step of the code we can then calculate the probability P tl that a planetesimal from bin i collides with any of the 
planetesimals in bin j during the time-step. If a random number r is less than / J , ; then the planetesimals will be collided. We 
assume perfect sticking between planetesimals and add the mass of the smaller planetesimal to the larger, removing afterwards 
the smaller planetesimal from the simulation. 

We tested the coagulation scheme against an analytical solution to the coagulation equation where the kernel is constant (62). 
This corresponds to setting the collision rate r l: j = Krij, where I\ is a constant and rij is the number density of planetesimals in 
bin j. The solution to the constant kernel test is that the number of bodies remaining at time t is 


n(t) 1 

no 1 + Knot/2 


(35) 


The number of remaining bodies versus time is shown in Figure |S9] for different values of the time-step. The time-step is set by 


d£ — CfltHcoll ; 


(36) 


where £ co n is the collision time-scale and Cdt is a constant time-step parameter. There is excellent agreement in Figure[S9|between 
the analytical expression and the results of the coagulation algorithm, for Cdt < 0.5. Passing this test implies that the Monte 
Carlo collision detection scheme is implemented correctly. 
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Figure S9: The number of remaining bodies versus time for the constant kernel test. There is an excellent agreement between 
the results of the coagulation algorithm and the analytical solution for a time-step parameter below cat = 0.5. The time-step 
parameter is the ratio of the time-step to the collision time-scale. 



t [yr] 


Figure S10: Damping of eccentricities and inclinations of 10,000 1-cm-sized planetesimals located between 30 and 35 AU, 
emulating a test problem defined in Morbidelli et al. (2009) (21). The total mass of the planetesimals is set to 10 Earth masses. 
The curves agree to within 15% with those shown in Figure 12 of the supplemental material of Morbidelli et al. (2009). 


A test of the collision rate calculation can not be done as simply, as the expressions used in Morbidelli et al. (2009) (21) are 
complex. We choose therefore to reproduce one of the figures in that paper, namely a plot of the damping of eccentricities and 
inclinations of planetesimals due to mutual inelastic collisions. The damping is not very relevant for the large planetesimals 
considered in this paper, but the exact shape of the damping provides an excellent comparison with the collision rate calculation 
of Morbidelli et al. (2009) (21). The test problem considers 10 Earth masses of 1-cm-sized “planetesimals” located between 30 
and 35 AU. Collisions are assumed to be inelastic and lead to energy dissipation; accretion and fragmentation are not included. 
Figure [ST0| shows the evolution of the eccentricity (initially 0.1) and the inclination (initially 0.05). There is agreement to within 
15% between this plot and supplemental Figure 12 of Morbidelli et al. (2009) (21). This shows that the collision rate algorithm 
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Figure SI 1: Cumulative size distribution after 3 Myr of coagulation within a population of planetesimals of initial diameters 100 
km (50 km in radius), as shown by the cross. The size distribution is much steeper than the current observed size distribution in 
the asteroid belt. The results presented here are agree to within a few percent with Figure 5 of Morbidelli et al. (2009) (27) and 
Figure 5 of Weidenschilling (2011) (57), except for D > 1000 km where we have a slight underproduction of embryos. Note 
that we set the column density of planetesimals equal to that of solids in the Minimum Mass Solar Nebula of Hayashi (1981), 
while Weidenschilling (2011) uses approximately twice that value. 


is likely implemented properly in the PAOPAP code. 

An additional test was performed with actual planetesimals. The evolution of a population of 50-km-radius planetesimals 
was presented in Figure 5 of Morbidelli et al. (2009) (27) and Figure 5 of Weidenschilling (2011) (57). We performed a 
similar simulation with our planetesimal accretion code. The resulting cumulative size distribution is shown in Figure |STT| The 
cumulative size distribution is very steep and terminates just before reaching embryo sizes. This steep cumulative size distribution 
agrees to within a few percent with the results of Morbidelli et al. (2009) (27) and Weidenschilling (2011) (57), except for 
a slightly smaller production of the very largest embryos in our simulations. Note that these authors include fragmentation 
of planetesimals in their calculations, while we ignore this effect. The agreement of the results for 79 > 20 km shows that 
planetesimal fragmentation plays little role for the coagulation of such large planetesimals. It is well-known that planetesimals 
larger than 100 km are strong and can survive high collision speeds (20). 
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Supplementary Text 


2.1 Planetesimal formation in the asteroid belt 

The particle sizes used in our streaming instability simulations correspond to approximately 25-cm-sized rocks at 2.5 AU in the 
Minimum Mass Solar Nebula (Stokes number St = 0.3). This is much larger than the typical diameters of chondrules. We 
have chosen to study the formation of planetesimals from such large particle sizes in order to have a clean convergence test in 
parameters employed in previous works (13,14). We suggest three ways to understand the actual formation of planetesimals in 
the asteroid belt: 

1) Chondrule aggregates. Chondrules are covered with rims of fine-grained matrix (4) and these rims can facilitate the 
formation of chondrule aggregates. If the turbulent speed of the gas is very low, with turbulent viscosity a ~ 10~ 6 , then 
collision speeds between rimmed chondrules are low enough to allow chondrules to stick together (77). The porous rims act 
as shock absorbers and facilitate the formation of decimeter-sized chondrule aggregates. These aggregates are large enough to 
trigger streaming instabilities and the direct formation of planetesimal seeds that can go on to accrete individual chondrules, as 
envisioned in the main paper. The radial drift speed of these aggregates is high, approximately 10 m/s, causing the aggregates to 
drift inwards within a few thousand years (63). The combination of growth and drift provides a means (a) to form planetesimals 
from such aggregates and (b) to subsequently flush these aggregates towards the inner solar system, to avoid feeding chondrule 
aggregates to the largest planetesimals, as this would trigger a very efficient run-away growth of asteroids and give a too steep 
size distribution. 

2) Large chondrules. Another possibility is that the first planetesimal seeds formed out of centimeter-sized macrochondrules 
(16). These could have drifted out of the asteroid belt subsequent to planetesimal formation and left the stage for accretion of 
their millimeter-sized counterparts. Carrera, Johansen, & Davies (45) investigate the ability for the streaming instability to form 
dense filaments in particles down to chondrule sizes. While chondrule-sized particles are very hard to concentrate, particles 
just a few times larger (Stokes number of 0.01) concentrate readily into filaments at a particle mass loading of a few times 
the nominal value in the solar protoplanetary disk (Zq ~ 0.01). Here we explore whether these conditions are susceptible to 
planetesimal formation, by including the self-gravity between centimeter-sized particles (Stokes number of 0.01, so 30 times 
smaller than in the simulations presented in Figure[3]of the main paper). We use 2-D shearing sheet simulations and set the mean 
particle density to 1, 3 and 10 times the gas density, respectively, to mimic the mid-plane conditions seen in Carrera, Johansen, 
& Davies (45). The results are shown in Figure |ST2| The simulations with particle density p p = ?>[>„ and p p = 10 p g form a 
single planetesimal after a few hundred years. These planetesimals have a characteristic radius of approximately 100 km. This 
shows that planetesimals can indeed form out of particles as small as 1 cm in size and that the resulting planetesimal sizes are 
similar to those seen in Figure |3]of the main paper. Thus we show that asteroids can form from centimeter-sized particles. The 
full exploration of the relevant parameters for formation of planetesimals from centimeter-sized particles is beyond the scope of 
this paper, but this should have a high priority in future research. 

3) Icy pebbles. It is possible that the ice line in the solar protoplanetary disk was much closer to the star in the earliest phases of 
protoplanetary disk evolution. Hence the asteroid belt could have been populated originally with icy pebbles of sizes comparable 
to what we used in the simulations (64). The first generation of icy asteroid seeds could subsequently have dried out after heating 
by 26 A1 or simply have been swamped in the subsequent accretion of dry chondrules. Carbonaceous chondrites display a wide 
range in their degree of aqueous alteration, consistent with the accretion and subsequent melting of icy pebbles (3). However, 
even some ordinary chondrites have been shown to have experienced water flows (65); and hence the formation of the seeds of 
the nowadays dry ordinary chondrites could have involved accretion of icy material. 

Irrespective of the details of asteroid formation, the important factors for our chondrule accretion model are (i) planetesimals 
form, (ii) the planetesimals are not too small to prevent significant chondrule accretion, and (iii) planetesimals are not so large 
that all chondrules are accreted on the largest objects in the population. The computer simulations of planetesimal formation 
by streaming instabilities inspire us to take asteroid seeds with characteristic radii of 100 km and with a fairly shallow size 
distribution that puts most of the mass in the largest bodies. 


2.2 Particle growth and chondrule precursors 

In our model we assume that the centimeter-to-decimeter-sized particles discussed in the previous section were only present 
during the earliest stages of planet formation when the planetesimals formed. Ice particles could only exist while the ice line was 
interior of the asteroid belt, while the formation of chondrule aggregates requires very weak turbulence of a ~ 10 -6 (77), which 
may have required an extensive dead zone in a young (and hence still massive) protoplanetary disk. In this section we discuss 
why the later stages of the protoplanetary disk were likely dominated by chondrule-sized particles. 
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Figure S12: Formation of planetesimals from centimeter-sized particles in a 2-D shearing sheet simulation with the same spatial 
extent as the streaming instability simulations presented in the main text. The mean particle density is p p = 1 in the left panel, 
p p = 3 in the middle panel and p p = 10 in the right panel (relative to the gas density). The mean density is set to mimic the 
effect of both stratification and filament formation by the streaming instability. The simulations with p p = 3 and p p = 10 form a 
single, large planetesimal of approximately 100 km in radius after a few hundred years of integration. This numerical experiment 
shows that planetesimals can form for particle sizes corresponding to centimeter-sized macrochondrules and that planetesimal 
sizes are similar to those forming out of 25-cm-sized particles presented in Figure[3]of the main text. 


The growth of dust particles in protoplanetary disks happens initially through sticking collisions that lead to the formation of 
fluffy dust aggregates; see the review paper on coagulation and planetesimal formation by Johansen et al. (7). These aggregates 
are compactified by mutual collisions as they reach approximately millimeter sizes. Compact dust aggregates have poor sticking 
properties and bounce off each other when they collide (66). This bouncing barrier is a major obstacle for dust growth (67), 
which could otherwise continue to much larger sizes (68,69). Particles can still grow by mass transfer in high-speed collisions; 
however the growth rates remain too low to compete against radial drift (70). Another possibility is that ice particles could be 
extremely fluffy and experience perfect sticking even at high speeds (77). This nevertheless requires aggregates to consist of very 
small monomers (0.1 /im), while matrix in chondrites is typically at least 10 times larger. 

Hence it appears that the growth of silicate dust grains in the asteroid formation region stops at around mm sizes. This bouncing 
barrier picture fits well with the sizes of chondrules found in chondrite meteorites. The actual formation process of chondrules 
is highly uncertain, but is believed to involve heating of dust aggregate precursors in shocks (72, 75), in current sheets (74) in 
the protoplanetary disk or in impacts between molten planetesimals (75). Our model does not include the presence of chondrule 
precursors. Such precursors could have a low fractal dimension and hence low friction time, which would prevent sedimentation 
and efficient incorporation into the asteroids and embryos that grow in the mid-plane. Nevertheless the inclusion of chondrule 
precursors would not change the conclusions of our paper, as they would simply lead to a slightly higher accretion rate for 
asteroids and embryos. 


2.3 Depletion of the asteroid belt 

The current mass of the asteroid belt is only approximately 5 x 10 -4 Earth masses. This contrasts with the extrapolation of the 
Minimum Mass Solar Nebula into the asteroid belt, which yields a particle column density of 4.3 g/cm 2 and hence a primordial 
mass of 0.86 Earth masses. The depletion of the asteroid belt has happened mainly through collisional grinding and pumping of 
the eccentricity of asteroids in resonances with Jupiter. Based on a number of arguments, including the extant number of asteroid 
families and the number of large craters on Vesta, Bottke et al. (20) concluded that the collisional activity in the asteroid belt, 
after excitation of the asteroids to their current eccentricities, can not have been much larger than the current collision rate over 
10 Gyr. Hence the asteroid belt must have been depleted to near its current mass very soon after asteroid formation. Morbidelli 
et al. (27) calculate that the asteroid belt could have lost only 1/3 of its current mass over a 10 Gyr equivalent of its present 
collisional activity. This low level of collisional grinding implies that the current size distribution of asteroids larger than 60 km 
in radius is primordial (20). 

The formation of embryos in the asteroid belt provides a means to excite the smaller asteroids to their current, high eccentric¬ 
ities (25). Embryos can also be responsible for the rapid clearing of the asteroid belt, since asteroids are gravitationally scattered 
by the embryos and hence move readily to resonances with Jupiter from which they are removed from the belt (76). Hence the 
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formation of planetary embryos in our chondrule accretion model is completely consistent with mechanisms for depletion of the 
asteroid belt. 

The main alternatives to the embryo model for asteroid belt depletion are (i) sweeping resonances (26) and (ii) the migration 
of Jupiter in the Grand Tack model (27). The migration of Jupiter and Saturn during the restructuring of the giant planets to their 
current orbits by scattering of the primordial Kuiper belt leads to a sweeping of resonant locations across the asteroid belt. This 
could have led to the loss of 90-95% of the mass of the asteroid belt during the migration period (26). The Grand Tack model 
provides a more violent means for Jupiter’s gravity to interfere with the asteroid belt, as Jupiter migrates through the asteroid belt 
twice; the change from inwards to outwards migration is caused by the formation of Saturn, which comes to share a gap with 
Jupiter in the protoplanetary disk. This injects C-type asteroids from beyond the ice line into the asteroid belt, which ends up 
strongly depleted; the later planetesimal-driven migration could inject further icy planetesimals into the main belt (77). 

Our simulations show that embryos form readily in the asteroid belt from the accretion of chondrules onto planetesimals. 
These embryos could have played an important role in sculpting the current orbits of the asteroids and depleting the asteroid 
belt. There are, as discussed above, mechanisms for the depletion of the asteroid belt that do not require the presence of embryos 
there. The reality could be that a combination of embedded embryos, sweeping resonances and Jupiter’s migration conspired to 
deplete the asteroid belt to its current mass. 
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